####setwd####

####packages####
options(scipen=999)

rm(list=ls())                    # Clears the memory.

# required packages
libraries<-c("xtable", "caseMatch","cowplot", "stargazer", "MatchIt",
             "foreign","raster", "readstata13",
             "dplyr",
             "gmodels", "lmtest", "maps","Matching","pastecs", "rgeos",
             "ggplot2", "readr", "rgenoud","rgdal",  "Rmisc", "stargazer", "tidyr")
lapply(libraries, library, character.only = TRUE)




# Department matching --------------------------------------------

# https://rdrr.io/cran/caseMatch/man/case.match.html
# Using caseMatch from  Nielsen, Richard A. 2016. “Case Selection via Matching.” Sociological Methods & Research 45(3):569–597.

dept_data<-read.csv("dept_data.csv")

#subset to matching variables
dept_data2<- dept_data %>% dplyr::select(dept, rebel, abj_dist, anti_i, north_1998, nightlights.2000.mean.x,
                                  africa_child_mortality.2000.mean.x)

# matching case selection
out <- case.match(data=dept_data2, id.var="dept",
                  distance="mahalanobis", case.N=2,
                  number.of.matches.to.return=15,
                  treatment.var="rebel", max.variance=TRUE)

# matching distances and cases
nielsen_match<-out$cases

# list of cases selected
nielsen_cases<-unique(c(unique(nielsen_match$`unit id 1`), unique(nielsen_match$`unit id 2`)))


# Table 1: Balance test of matched pairs ----------------------------------

# subset the dept_data to the matched sample
nielsen_match_data<- dept_data %>% filter(dept %in% nielsen_cases)

match_vars<-c("abj_dist", "anti_i", "nightlights.2000.mean.x",
              "africa_child_mortality.2000.mean.x", "north_1998")

match_vars2<-c("Distance to Abidjan, km", "Anti-Incumbent, 2001", "Night Lights, 2000", "Child Mortality, 2000",
               "Northern Ethnic Groups, 1998")

# create top panel of Table 1 with the matching variables balance

table_1a<-data.frame(match_vars2)

for (i in 1:length(match_vars)){
  table_1a[i,2]<-lapply(nielsen_match_data[match_vars[i]], function(x) {t.test(x[nielsen_match_data$rebel==1], x[nielsen_match_data$rebel==0])$estimate[1]})
  table_1a[i,3]<-lapply(nielsen_match_data[match_vars[i]], function(x) {t.test(x[nielsen_match_data$rebel==1], x[nielsen_match_data$rebel==0])$estimate[2]})
  table_1a[i,4]<-lapply(nielsen_match_data[match_vars[i]], function(x) {t.test(x[nielsen_match_data$rebel==1], x[nielsen_match_data$rebel==0])$p.value})
}

table_1a$N<-nrow(nielsen_match_data)

names(table_1a)<-c("Variable", "Mean Rebel", "Mean Govt", "P-value", "N")
table_1a$Diff<-table_1a$`Mean Rebel`- table_1a$`Mean Govt`
table_1a<- table_1a %>% dplyr::select(Variable,  `Mean Rebel`, `Mean Govt`, Diff, `P-value`, N) %>% mutate_if(is.numeric, round, 2)

print(xtable(table_1a), include.rownames=F, file="table_1a.tex")

# create panel 2 of Table 1 with the other pre-war variables
prewar_vars<-c("akan_1998",
               "krou_1998",
               "smande_1998",
               "muslim_98",
               "X1998_pop",
               "land_inegal_1998",
               "health_pop_b42002",
               "avg_rain",
               "cacao_production",
               "distance_to_gold_v12.none.mean.x",
               "prewar_battles", "min_liberia_border")

prewar_vars2<-c("Akan, 1998",
                "Krou, 1998",
                "South Mandé, 1998",
                "Muslim, 1998",
                "1998 Population",
                "Land inequality, gini coefficient",
                "Health centers per 1000 inhabitants, 2001",
                "Average 10-year rainfall (1992-2002), in",
                "Cacao production, kg",
                "Distance to lootable gold, km",
                "Pre-war battles at least one in dept, 1992-2001", "Distance to Liberian border")


table_1b<-data.frame(prewar_vars2)

for (i in 1:length(prewar_vars)){
  table_1b[i,2]<-lapply(nielsen_match_data[prewar_vars[i]], function(x) {t.test(x[nielsen_match_data$rebel==1], x[nielsen_match_data$rebel==0])$estimate[1]})
  table_1b[i,3]<-lapply(nielsen_match_data[prewar_vars[i]], function(x) {t.test(x[nielsen_match_data$rebel==1], x[nielsen_match_data$rebel==0])$estimate[2]})
  table_1b[i,4]<-lapply(nielsen_match_data[prewar_vars[i]], function(x) {t.test(x[nielsen_match_data$rebel==1], x[nielsen_match_data$rebel==0])$p.value})
}

table_1b$N<-nrow(nielsen_match_data)

names(table_1b)<-c("Variable", "Mean Rebel", "Mean Govt", "P-value", "N")
table_1b$Diff<-table_1b$`Mean Rebel`- table_1b$`Mean Govt`
table_1b<- table_1b %>% dplyr::select(Variable,  `Mean Rebel`, `Mean Govt`, Diff, `P-value`, N)  %>% mutate_if(is.numeric, round, 2)

print(xtable(table_1b), include.rownames=F, file="table_1b.tex")


# Figure 1: Selected cases in Côte d’Ivoire -------------------------------

# read in csv from shape files
civ.df<-read_csv("civ_df.csv")
names(civ.df)
dept_centroids<-read_csv("civ_dept_centroids_FINAL.csv")
table(civ.df$cacao)
#read in zdc csv
zdc_upper<-read_csv("zdc_upper.csv")
zdc_lower<-read_csv("zdc_lower.csv")

# read in neighboring country data to make map

liberia = readRDS("GADM_2.8_LBR_adm0.rds", refhook = NULL)
guinea = readRDS("GADM_2.8_GIN_adm0.rds", refhook = NULL)
mali = readRDS("GADM_2.8_MLI_adm0.rds", refhook = NULL)
ghana =readRDS("GADM_2.8_GHA_adm0.rds", refhook = NULL)
burkina = readRDS("GADM_2.8_BFA_adm0.rds", refhook = NULL)


#this is just so there isn't overlap in the name of the deparment and the label for the labs
dept_centroids_cases<-dept_centroids %>% filter(dept_name %in% nielsen_cases) %>% filter(dept_name!="San Pedro") %>% filter(dept_name!="Man") %>%
  filter(dept_name!="Daloa")


# getting shape files for matched cases
case_map<-civ.df %>% mutate(cases=ifelse(civ.df$NAME_3 %in% nielsen_cases, 1, 0)) %>% filter(cases==1)


# lab locations
lab<- dept_centroids %>% mutate(lab=ifelse(dept_name=="San Pedro" |
                             dept_name=="Daloa"| dept_name=="Man"|
                             dept_name=="Bouake", 1, 0)) %>% filter(lab==1)

# capitals

capitals<- dept_centroids %>% mutate(capitals=ifelse(dept_name=="Abidjan" | dept_name=="Yamoussoukro", 1, 0)) %>%
  filter(capitals==1)



# Map!
ggplot(civ.df)+
  aes(x=long, y=lat) +
  geom_polygon(aes(group=group, fill=factor(civ.df$rebel)), color="grey")+
  scale_fill_manual(values=c("white", "darkgray"),
                    labels=c("Government controlled ", "Rebel controlled"))+
  geom_text(data=dept_centroids_cases, aes(label = dept_name, x = long, y = lat), color="black", size=3)+
  geom_point(data=capitals,aes(x=long,y=lat), fill="black", shape=17, size=2)+
  geom_text(data=capitals, aes(label = dept_name, x = long, y = lat),
            color="black", size=3.5,nudge_y = -0.04,  nudge_x = .3, family="serif")+
  geom_path(data=case_map, aes(x=long, y=lat, group=group), color="black", size=1)+
  geom_line(data=zdc_lower, col="black",linetype = 2, size=1)+
  geom_line(data=zdc_upper, col="black",linetype = 2, size=1)+
  geom_polygon(data=liberia,fill=NA,color="black")+
  geom_polygon(data=guinea,fill=NA,color="black")+
  geom_polygon(data=mali,fill=NA,color="black")+
  geom_polygon(data=ghana,fill=NA,color="black")+
  geom_polygon(data=burkina,fill=NA,color="black")+
  coord_map(xlim=c(-9.599302,-1.494897),ylim=c(3.361807,11.73664))+
  theme_bw()+theme(panel.grid=element_blank(), axis.text.x=element_blank(),
                   axis.text.y=element_blank(), axis.ticks.x=element_blank(),
                   axis.ticks.y=element_blank(),
                   axis.title.x=element_blank(), axis.title.y=element_blank())+
  annotate("text", x = -9, y = 6, label = "Liberia")+
  annotate("text", x = -9, y = 9, label = "Guinea")+
  annotate("text", x = -7, y = 11, label = "Mali")+
  annotate("text", x = -4, y = 11, label = "Burkina Faso")+
  annotate("text", x = -2, y = 7, label = "Ghana")+
  guides(fill="none")+
  geom_label(data=lab, aes(label = dept_name, x = long, y = lat), color="black",
             size=4, nudge_y = -0.1, nudge_x=0.05, family="serif")


ggsave(filename="Fig1.pdf", device = cairo_pdf, width = 15, height=8.5, units="in", dpi = 600)

# Table A1: balance on matching covariates, pre-treatment -----------------


Table_A1<-data.frame(match_vars2)

for (i in 1:length(match_vars)){
  Table_A1[i,2]<-lapply(dept_data[match_vars[i]], function(x) {t.test(x[dept_data$rebel==1], x[dept_data$rebel==0])$estimate[1]})
  Table_A1[i,3]<-lapply(dept_data[match_vars[i]], function(x) {t.test(x[dept_data$rebel==1], x[dept_data$rebel==0])$estimate[2]})
  Table_A1[i,4]<-lapply(dept_data[match_vars[i]], function(x) {t.test(x[dept_data$rebel==1], x[dept_data$rebel==0])$p.value})
}

Table_A1$N<-nrow(dept_data)

names(Table_A1)<-c("Variable", "Mean Rebel", "Mean Govt", "P-value", "N")
Table_A1$Diff<-Table_A1$`Mean Rebel`- Table_A1$`Mean Govt`
Table_A1<- Table_A1 %>% dplyr::select(Variable,  `Mean Rebel`, `Mean Govt`, Diff, `P-value`, N) %>% mutate_if(is.numeric, round, 2)

print(xtable(caption="Balance on matching covariates, pre-treatment", Table_A1), include.rownames=F, file="Table_A1.tex")


# Table A2: Cases using Mahalanobis distance matching -------------------------------------------------------------------------

nielsen_match$pair<-rownames(nielsen_match)
Table_A2<-nielsen_match %>% pivot_longer(`unit id 1`:`unit id 2`) %>% dplyr::select(pair, value, distances)
names(Table_A2)<-c("Pair", "dept", "d")

nielsen_match_data_tblA2<- nielsen_match_data %>% dplyr::select(dept, rebel, abj_dist,
                                                                nightlights.2000.mean.x,
                                                                africa_child_mortality.2000.mean.x,
                                                                north_1998,
                                                                anti_i)

Table_A2<-merge(Table_A2, nielsen_match_data_tblA2, by="dept")
Table_A2<- Table_A2 %>% dplyr::select(Pair, dept, rebel, abj_dist, nightlights.2000.mean.x, africa_child_mortality.2000.mean.x, north_1998, anti_i, d) %>%
  arrange(d) %>% mutate_if(is.numeric, round, 2)
names(Table_A2)<-c("Pair", "Dept", "Rebel", "Dist. Abidjan", "Night Light 2000", "Child Mort. 2000", "Northern 1998", "Anti-Incumbent 2001", "d")

# To make this table in the Appendix, I hand removed the duplicates of "Pair" and "d" and added horizontal lines at each pair break.
print(xtable(Table_A2), include.rownames=F, file="Table_A2.tex")


# Table A3: Balance test of matched pairs: post treatment variables -------

# creating top panel for A3 for intra-war variables balance
intra_vars<- c("vac.cw",
"battles.cw",
"projects20022008",
"USD20022008",
"huma_2006_count",
"vote_2010_r2_gbagbo",
"vac.ev",
"battles.ev")


intra_vars2<-c("Violence against civilians, 2002-2008","Battles, 2002-2008",
               "Count of aid projects, 2002-2005",
               "Constant USD of aid projects, 2002-2005",
               "Number of int. aid orgs, 2006",
               "R2 Incumbent vote share, 2010",
               "Violence against civilians, 2010-2011",
               "Battles, 2010-2011")


table_A3a<-data.frame(intra_vars2)

for (i in 1:length(intra_vars)){
  table_A3a[i,2]<-lapply(nielsen_match_data[intra_vars[i]], function(x) {t.test(x[nielsen_match_data$rebel==1], x[nielsen_match_data$rebel==0])$estimate[1]})
  table_A3a[i,3]<-lapply(nielsen_match_data[intra_vars[i]], function(x) {t.test(x[nielsen_match_data$rebel==1], x[nielsen_match_data$rebel==0])$estimate[2]})
  table_A3a[i,4]<-lapply(nielsen_match_data[intra_vars[i]], function(x) {t.test(x[nielsen_match_data$rebel==1], x[nielsen_match_data$rebel==0])$p.value})
}

table_A3a$N<-nrow(nielsen_match_data)

names(table_A3a)<-c("Variable", "Mean Rebel", "Mean Govt", "P-value", "N")
table_A3a$Diff<-table_A3a$`Mean Rebel`- table_A3a$`Mean Govt`
table_A3a<- table_A3a %>% dplyr::select(Variable,  `Mean Rebel`, `Mean Govt`, Diff, `P-value`, N) %>% mutate_if(is.numeric, round, 2)


print(xtable(table_A3a), include.rownames=F, file="table_A3a.tex")

# creating top panel for A3 for post-war variables balance
post_war_vars<-c("nightlights.2016.mean",
"X2014_pop",
"akan_2014",
"krou_2014",
"SMande_2014",
"northern_2014")

post_war_vars2<-c("Night lights, 2016",
                  "2014 population",
                  "Akan, 2014",
                  "Krou, 2014",
                  "South Mandé, 2014",
                  "Northern ethnic groups, 2014")

table_A3b<-data.frame(post_war_vars2)

for (i in 1:length(post_war_vars)){
  table_A3b[i,2]<-lapply(nielsen_match_data[post_war_vars[i]], function(x) {t.test(x[nielsen_match_data$rebel==1], x[nielsen_match_data$rebel==0])$estimate[1]})
  table_A3b[i,3]<-lapply(nielsen_match_data[post_war_vars[i]], function(x) {t.test(x[nielsen_match_data$rebel==1], x[nielsen_match_data$rebel==0])$estimate[2]})
  table_A3b[i,4]<-lapply(nielsen_match_data[post_war_vars[i]], function(x) {t.test(x[nielsen_match_data$rebel==1], x[nielsen_match_data$rebel==0])$p.value})
}

table_A3b$N<-nrow(nielsen_match_data)

names(table_A3b)<-c("Variable", "Mean Rebel", "Mean Govt", "P-value", "N")
table_A3b$Diff<-table_A3b$`Mean Rebel`- table_A3b$`Mean Govt`
table_A3b<- table_A3b %>% dplyr::select(Variable,  `Mean Rebel`, `Mean Govt`, Diff, `P-value`, N) %>% mutate_if(is.numeric, round, 2)


print(xtable(table_A3b), include.rownames=F, file="table_A3b.tex")


# Figure A1: The location of the labs and the departments --------


ggplot(civ.df)+
  aes(x=long, y=lat) +
  geom_polygon(aes(group=group, fill=civ.df$lab), color='gray')+
  coord_map("mercator")+
  theme_classic()+
  scale_fill_brewer(name= "Lab locations, cases", palette="Dark2", na.translate = F)+
  geom_path(data=case_map, aes(x=long, y=lat, group=group), color="black", size=1)+
  geom_label(data=lab, aes(label = dept_name, x = long, y = lat), color="black",
             size=5, nudge_y = -0.05, nudge_x=-0.05, family="serif")+
  geom_line(data=zdc_lower, col="black",linetype = 2, size=1)+
  geom_line(data=zdc_upper, col="black",linetype = 2, size=1)+
  theme(legend.text = element_text(size=20), legend.title = element_text(size=20))

ggsave(filename="FigA1.pdf", device = cairo_pdf, width = 15, height=8.5, units="in", dpi = 600)


# LAB-IN-THE-FIELD ANALYSIS --------------------------------------------------------------------------
labs<-read.csv("labs_df.csv", stringsAsFactors=FALSE, as.is=T)

# Table 2: Summary Statistics for lab-in-the-field sample -----------------------------------------------------------------


dem<-c("rebel", "SEX", "over35", "higher_edu", "muslim", "northern", "support_incumbent", "INCOME2", "remained",
       "left",
       "FOUNDER", "prior_ngo", "created_prewar", "created_war","int_fund", "natl", "total_domains",
       "DOMAIN_HUMA", "ORG_MISSION_INFO_2")


description<-c("CSO leader from rebel depts" , "Female", "Age above 35", "Some secondary education",
               "Muslim","Northern ethnic group (North Mandé or Voltaïque)",
               "Party id is president's party",
               "Income covers expenses", "Did not migrate during war",
               "Migrated during war",
               "Founder", "Worked in NGOs prior to current NGO",
               "Created before war", "Created during war","Internationally funded", "Natl headquarters",
               "Total number of domains","Humanitarian","Info provision")

Table2<-t(stat.desc(labs[dem]))

Table2<-as.data.frame(Table2[,c("mean", "min", "max", "std.dev", "nbr.val")])

row.names(Table2)<-description
colnames(Table2)<-c("mean", "min", "max", "std.dev", "n")
Table2<- Table2 %>% mutate_if(is.numeric, round, 2)

print(xtable(Table2), file="Table2.tex")

# Table 3: Outcome summary statistics -------------------------------------

outcomes<-c("self_ind", "coethnic", "noncoethnic", "copartisan", "noncopartisan", "outgroup")
outcomes2<-c("Self", "Coethnic", "Noncoethnic", "Copartisan","Noncopartisan", "Outgroup")

Table3<-t(stat.desc(labs[outcomes]))

Table3<-as.data.frame(Table3[,c("mean", "min", "max", "std.dev", "nbr.val")])

row.names(Table3)<-outcomes2
colnames(Table3)<-c("mean", "min", "max", "std.dev", "n")

Table3<- Table3 %>% dplyr::select(mean, min, max, std.dev, n) %>% mutate_if(is.numeric, round, 2)

print(xtable(Table3), file="Table3.tex")


# Figure 3: Altruism: Lower redistribution in rebel zones -----------------

self_plot<-summarySE(labs, measurevar="self_ind_pct", groupvars = c("rebel"), na.rm=T)


Fig3a<-ggplot()+
  geom_point(data=labs, aes(x=factor(rebel), y=100*self_ind_pct), alpha = .1)+
  #scale_x_discrete(labels=c("Govt", "Rebel"))+
  geom_bar(data=self_plot, aes(x=factor(rebel), y=100*self_ind_pct, fill=factor(rebel)), alpha=.5, stat="identity", width=.5, color="black")+
  geom_errorbar(data=self_plot, aes(x=factor(rebel), y=100*self_ind_pct, ymin=100*(self_ind_pct-ci), ymax=100*(self_ind_pct+ci)),
                width = 0.1, position = position_dodge(width = 1/2)) +
  theme_minimal()+
  geom_segment(aes(x=1, y=100, xend=1.1, yend=103), colour="darkgray")+
  geom_segment(aes(x=1.1, y=103, xend=1.9, yend=103), colour="darkgray")+
  geom_segment(aes(x=1.9, y=103, xend=2, yend=100), colour="darkgray")+
  annotate("text", x=1.5, y=105, label="11%p (p<0.001)", size=7)+
  scale_x_discrete(labels=c("Leader from \nformer-Govt", "Leader from \nformer-Rebel"))+
  scale_fill_manual(values=c("white", "darkgray"))+
  theme(legend.position ="none", axis.text  = element_text(size=20),
        plot.title = element_text(size = 20, hjust = 0.5, face = "bold"),
        panel.border = element_rect(colour = "black", fill=NA))+
  labs(y = "", x="")+
  ggtitle("Allocate to self (%)")

out_plot_total<-summarySE(labs, measurevar="outgroup_pct", groupvars = c("rebel"), na.rm=T)


Fig3b<-ggplot()+
  geom_point(data=labs, aes(x=factor(rebel), y=100*outgroup_pct), alpha = .1)+
  #scale_x_discrete(labels=c("Govt", "Rebel"))+
  geom_bar(data=out_plot_total, aes(x=factor(rebel), y=100*outgroup_pct, fill=factor(rebel)), alpha=.5, stat="identity", width=.5, color="black")+
  geom_errorbar(data=out_plot_total, aes(x=factor(rebel), y=100*outgroup_pct, ymin=100*(outgroup_pct-ci), ymax=100*(outgroup_pct+ci)),
                width = 0.1, position = position_dodge(width = 1/2)) +
  #geom_point(data=out_plot_total2, aes(x=rebel, y=100*outgroup_pct, fill=rebel), stat="identity", size=4, shape=24, stroke=2)+
  scale_fill_manual(values=c("white", "darkgray"))+
  theme_minimal()+
  scale_x_discrete(labels=c("Leader from \nformer-Govt", "Leader from \nformer-Rebel"))+
  geom_segment(aes(x=1, y=100, xend=1.1, yend=103), colour="darkgray")+
  geom_segment(aes(x=1.1, y=103, xend=1.9, yend=103), colour="darkgray")+
  geom_segment(aes(x=1.9, y=103, xend=2, yend=100), colour="darkgray")+
  annotate("text", x=1.5, y=105, label="7%p (p<0.001)", size=7)+
  theme(legend.position ="none", axis.text  = element_text(size=20),
        plot.title = element_text(size = 20, hjust = 0.5, face = "bold"),
        axis.text.y=element_blank(),
        panel.border = element_rect(colour = "black", fill=NA))+
  labs(y = "", x="")+
  scale_y_continuous(limits=c(0,105))+
  ggtitle("Allocate to outgroup (%)")


plot_grid(Fig3a, Fig3b)
ggsave(filename="Fig3.pdf", device = cairo_pdf, width = 15, height=8.5, units="in", dpi = 600)

# Table 4: Egocentrism and Discrimination among Civil Society Lead --------


base_self.lm<-lm(self_ind_pct~rebel, data=labs)
base_self.lm2<-coeftest(base_self.lm, vcov =  sandwich::vcovCL(base_self.lm, labs$DEPT[is.na(labs$rebel) == FALSE]), df = 20)

self_ind.lm<-lm(self_ind_pct~rebel+SEX+higher_edu+northern, data=labs)
self_ind.lm2<-coeftest(self_ind.lm, vcov =  sandwich::vcovCL(self_ind.lm, labs$DEPT[is.na(labs$rebel) == FALSE  &
                                                                                     is.na(labs$SEX) == FALSE &
                                                                                     is.na(labs$self_ind_pct)==F &
                                                                                     is.na(labs$higher_edu)==F &
                                                                                     is.na(labs$northern)==F]), df = 20)
org_ctrls<-lm(self_ind_pct~rebel+int_fund+natl+total_domains+created_war+created_postwar+DOMAIN_HUMA+ORG_MISSION_INFO_2, data=labs)

org_ctrls.self<-coeftest(org_ctrls, vcov =  sandwich::
                           vcovCL(org_ctrls, labs$DEPT.x[is.na(labs$rebel) == FALSE  &
                                                             is.na(labs$int_fund) == FALSE &
                                                             is.na(labs$self_ind_pct)==F &
                                                             is.na(labs$natl)==F &
                                                             is.na(labs$total_domains)==F &
                                                             is.na(labs$created_postwar)==F &
                                                             is.na(labs$DOMAIN_HUMA)==F &

                                                             is.na(labs$ORG_MISSION_INFO_2)==F &
                                                             is.na(labs$created_war)==F]), df = 20)


base_out.lm<-lm(outgroup_pct~rebel, data=labs)
base_out.lm2<-coeftest(base_out.lm, vcov =  sandwich::vcovCL(base_out.lm, labs$DEPT[is.na(labs$rebel) == FALSE]), df = 20)

out.lm<-lm(outgroup_pct~rebel+SEX+higher_edu+northern, data=labs)
out.lm2<-coeftest(out.lm, vcov =  sandwich::vcovCL(out.lm, labs$DEPT[is.na(labs$rebel) == FALSE  &
                                                                      is.na(labs$SEX) == FALSE &
                                                                      is.na(labs$outgroup_pct)==F &
                                                                      is.na(labs$higher_edu)==F &
                                                                      is.na(labs$northern)==F]), df = 20)
org_ctrls_out<-lm(outgroup_pct~rebel+int_fund+natl+total_domains+created_war+created_postwar+DOMAIN_HUMA+ORG_MISSION_INFO_2, data=labs)
org_ctrls_out.out<-coeftest(org_ctrls_out, vcov =  sandwich::
                              vcovCL(org_ctrls_out, labs$DEPT.x[is.na(labs$rebel) == FALSE  &
                                                                    is.na(labs$int_fund) == FALSE &
                                                                    is.na(labs$outgroup_pct)==F &
                                                                    is.na(labs$natl)==F &
                                                                    is.na(labs$total_domains)==F &
                                                                    is.na(labs$created_postwar)==F &
                                                                    is.na(labs$DOMAIN_HUMA)==F &

                                                                    is.na(labs$ORG_MISSION_INFO_2)==F &
                                                                    is.na(labs$created_war)==F]), df = 20)



full_reg<-lm(self_ind_pct~rebel+int_fund+natl+total_domains+created_war+created_postwar+DOMAIN_HUMA+ORG_MISSION_INFO_2+SEX+higher_edu+northern, data=labs)
full_reg_self<-coeftest(full_reg, vcov =  sandwich::
                     vcovCL(full_reg, labs$DEPT.x[is.na(labs$rebel) == FALSE  &
                                                      is.na(labs$int_fund) == FALSE &
                                                      is.na(labs$self_ind_pct)==F &
                                                      is.na(labs$natl)==F &
                                                      is.na(labs$total_domains)==F &
                                                      is.na(labs$created_postwar)==F &
                                                      is.na(labs$DOMAIN_HUMA)==F &
                                                      is.na(labs$SEX) == FALSE &
                                                      is.na(labs$higher_edu)==F &
                                                      is.na(labs$northern)==F &
                                                      is.na(labs$ORG_MISSION_INFO_2)==F &
                                                      is.na(labs$created_war)==F]), df = 20)

full_reg2<-lm(outgroup_pct~rebel+int_fund+natl+total_domains+created_war+created_postwar+DOMAIN_HUMA+ORG_MISSION_INFO_2+SEX+higher_edu+northern, data=labs)
full_reg_out<-coeftest(full_reg2, vcov =  sandwich::
                      vcovCL(full_reg2, labs$DEPT.x[is.na(labs$rebel) == FALSE  &
                                                        is.na(labs$int_fund) == FALSE &
                                                        is.na(labs$outgroup_pct)==F &
                                                        is.na(labs$natl)==F &
                                                        is.na(labs$total_domains)==F &
                                                        is.na(labs$created_postwar)==F &
                                                        is.na(labs$DOMAIN_HUMA)==F &
                                                        is.na(labs$SEX) == FALSE &
                                                        is.na(labs$higher_edu)==F &
                                                        is.na(labs$northern)==F &
                                                        is.na(labs$ORG_MISSION_INFO_2)==F &
                                                        is.na(labs$created_war)==F]), df = 20)





require(stargazer)
stargazer(base_self.lm2, self_ind.lm2, org_ctrls.self, full_reg_self,
          base_out.lm2, out.lm2, org_ctrls_out.out,full_reg_out, style = "apsr",
          type="text",
          column.labels   = c("Percent allocated to self", "Percent allocated to outgroups"),
          column.separate = c(4, 4),
          add.lines = list(c("Observations", paste(nobs(base_self.lm)),paste(nobs(self_ind.lm)),
                             paste(nobs(org_ctrls)),paste(nobs(full_reg2)),
                             paste(nobs(base_out.lm)),paste(nobs(out.lm)),
                             paste(nobs(org_ctrls_out)),paste(nobs(full_reg2))),
                           c("Individual controls", "No", "Yes","No", "Yes","No", "Yes", "No", "Yes"),
                           c("Organization controls", "No", "No","Yes", "Yes","No", "No", "Yes", "Yes")),
          keep = c("rebel"),
          covariate.labels = "Civil society leader from rebel depts",
          out = "Table4.tex")

# Table 5: Individual experience with uncertainty -------------------------

mech<-c("FEAR_20022008_DANGER", "FEAR_20022008_RULEOFLAW",
       "FEAR_20022008_UNKNOWN",
       "FEAR_20022008_REBELS",
       "FEAR_20022008_GOVT",
        "feared_discrimination")

mech2<-c("I felt in danger", "I felt rule of law persisted","I didn't know who ruled", "I feared the rebels",
         "I feared the govt", "I faced discrimination (ethnic, religious, or partisanship)")
Table5<-data.frame(mech2)


for (i in 1:length(mech)){
  Table5[i,2]<-lapply(labs[mech[i]], function(x) {t.test(x[labs$rebel==1], x[labs$rebel==0])$estimate[1]})
  Table5[i,3]<-lapply(labs[mech[i]], function(x) {t.test(x[labs$rebel==1], x[labs$rebel==0])$estimate[2]})
  Table5[i,4]<-lapply(labs[mech[i]], function(x) {t.test(x[labs$rebel==1], x[labs$rebel==0])$p.value})
  Table5[i,5]<-lapply(labs[mech[i]], function(x) {min(x, na.rm=T)})
  Table5[i,6]<-lapply(labs[mech[i]], function(x) {max(x, na.rm=T)})
  Table5[i,7]<-lapply(labs[mech[i]], function(x) {sum(!is.na(x))})

}


names(Table5)<-c("Variable", "Mean Rebel", "Mean Govt", "P-value", "min","max", "N")
Table5$Diff<-Table5$`Mean Rebel`- Table5$`Mean Govt`
Table5<- Table5 %>% dplyr::select(Variable,  `Mean Rebel`, `Mean Govt`, Diff, `P-value`, min, max, N)%>% mutate_if(is.numeric, round, 2)

print(xtable(Table5), include.rownames=F, file="Table5.tex")


# Table A4: balance on individual-level variables -------------------------

#here include those in table 2 + additional vars

ind_balance<-c("SEX", "over35", "higher_edu", "muslim", "northern", "support_incumbent", "INCOME2",
               "remained", "arrived",
               "left", "EGOTISM2")

ind_balance2<-c("Female", "Over 35years", "Some secondary education", "Muslim", "Northerner",
                "Party ID is president's party", "Income covers expenses","Did not migrate during war",
                "Arrived during war", "Migrated during war", "Only help if I benefit")

TableA4<-data.frame(ind_balance2)

for (i in 1:length(ind_balance)){
  TableA4[i,2]<-lapply(labs[ind_balance[i]], function(x) {t.test(x[labs$rebel==1], x[labs$rebel==0])$estimate[1]})
  TableA4[i,3]<-lapply(labs[ind_balance[i]], function(x) {t.test(x[labs$rebel==1], x[labs$rebel==0])$estimate[2]})
  TableA4[i,4]<-lapply(labs[ind_balance[i]], function(x) {t.test(x[labs$rebel==1], x[labs$rebel==0])$p.value})
  TableA4[i,5]<-lapply(labs[ind_balance[i]], function(x) {sum(!is.na(x))})
  TableA4[i,6]<-lapply(labs[ind_balance[i]], function(x) {min(x, na.rm=T)})
  TableA4[i,7]<-lapply(labs[ind_balance[i]], function(x) {max(x, na.rm=T)})

}

names(TableA4)<-c("Variable", "Mean Rebel", "Mean Govt", "P-value", "N", "min","max")
TableA4$Diff<-TableA4$`Mean Rebel`- TableA4$`Mean Govt`
TableA4<- TableA4 %>% dplyr::select(Variable,  `Mean Rebel`, `Mean Govt`, Diff, `P-value`,min, max, N) %>% mutate_if(is.numeric, round, 2)

print(xtable(TableA4), include.rownames=F, file="TableA4.tex")

# Table A5: Balance on NGO-level variables --------------------------------

#here include those in table 2 + all domains + partnerships

org_balance<-c("FOUNDER", "prior_ngo", "created_prewar", "created_war","int_fund", "natl",
               "EXECUTED_c", "MEMBER_COUNT_c", "total_domains",
               "ORG_MISSION_INFO_2", "ORG_MISSION_PRIVATE","ORG_MISSION_PUBLIC" ,
               "local_part", "min_part", "rural")

org_balance2<-c("Founder", "Worked in NGOs prior to current NGO",
                "Created before war", "Created during war","Internationally funded",
                "Natl headquarters",
                "Projects (count)",
                "Number of members", "Total domains",
                "Mission: Information provision",
                "Mission: private goods provision",
                "Mission: Public goods provision",
                "Local government partnerships",
                "National government partnerships",
                "Rural projects")

TableA5<-data.frame(org_balance2)

for (i in 1:length(org_balance)){
  TableA5[i,2]<-lapply(labs[org_balance[i]], function(x) {t.test(x[labs$rebel==1], x[labs$rebel==0])$estimate[1]})
  TableA5[i,3]<-lapply(labs[org_balance[i]], function(x) {t.test(x[labs$rebel==1], x[labs$rebel==0])$estimate[2]})
  TableA5[i,4]<-lapply(labs[org_balance[i]], function(x) {t.test(x[labs$rebel==1], x[labs$rebel==0])$p.value})
  TableA5[i,5]<-lapply(labs[org_balance[i]], function(x) {sum(!is.na(x))})
  TableA5[i,6]<-lapply(labs[org_balance[i]], function(x) {min(x, na.rm=T)})
  TableA5[i,7]<-lapply(labs[org_balance[i]], function(x) {max(x, na.rm=T)})

}


names(TableA5)<-c("Variable", "Mean Rebel", "Mean Govt", "P-value", "N", "min","max")
TableA5$Diff<-TableA5$`Mean Rebel`- TableA5$`Mean Govt`
TableA5<- TableA5 %>% dplyr::select(Variable,  `Mean Rebel`, `Mean Govt`, Diff, `P-value`,  min, max,N) %>% mutate_if(is.numeric, round, 2)

print(xtable(TableA5), include.rownames=F, file="TableA5.tex")


# Table A6: Domains balance -------------------------------------------

### show t-tests for these domains
domains.df<- labs %>% filter(labs$ResponseId!="NOT IN DATASET")

domains_top<-c("DOMAIN_DH",
"DOMAIN_HUMA",
"DOMAIN_COHESION",
"DOMAIN_CHILD",
"DOMAIN_DEV",
"DOMAIN_FEMME",
"DOMAIN_HEALTH")

domains_top2<-c("Human rights","Humanitarian","Social Cohesion","Children","Development","Women","Health")

TableA6<-data.frame(domains_top2)

for (i in 1:length(domains_top)){
  TableA6[i,2]<-lapply(domains.df[domains_top[i]], function(x) {t.test(x[domains.df$rebel==1], x[domains.df$rebel==0])$estimate[1]})
  TableA6[i,3]<-lapply(domains.df[domains_top[i]], function(x) {t.test(x[domains.df$rebel==1], x[domains.df$rebel==0])$estimate[2]})
  TableA6[i,4]<-lapply(domains.df[domains_top[i]], function(x) {t.test(x[domains.df$rebel==1], x[domains.df$rebel==0])$p.value})
  TableA6[i,5]<-lapply(domains.df[domains_top[i]], function(x) {sum(!is.na(x))})
  TableA6[i,6]<-lapply(domains.df[domains_top[i]], function(x) {min(x, na.rm=T)})
  TableA6[i,7]<-lapply(domains.df[domains_top[i]], function(x) {max(x, na.rm=T)})

}


names(TableA6)<-c("Variable", "Mean Rebel", "Mean Govt", "P-value", "N", "min","max")

#selecting domains for which there is not balance
domain_imbal<-TableA6$Variable[TableA6$`P-value`<=.05][-1]


TableA6$Diff<-TableA6$`Mean Rebel`- TableA6$`Mean Govt`
TableA6<- TableA6 %>% dplyr::select(Variable,  `Mean Rebel`, `Mean Govt`, Diff, `P-value`, min, max, N) %>% mutate_if(is.numeric, round, 2)

print(xtable(TableA6), include.rownames=F, file="TableA6.tex")

# Now creating table for the outcomes for imbalanced domains
TableA7a<-data.frame(domain_imbal)

# top panel of A7
domains_top_sig<-c("DOMAIN_HUMA","DOMAIN_CHILD","DOMAIN_FEMME")
for (i in 1:length(domains_top_sig)){
  tryCatch({
    TableA7a[i,2]<-lapply(domains.df[domains_top_sig[i]], function(x) {t.test(domains.df$self_ind_pct[x==1 & domains.df$rebel==1], domains.df$self_ind_pct[x==1 & domains.df$rebel==0])$estimate[1]})
    TableA7a[i,3]<-lapply(domains.df[domains_top_sig[i]], function(x) {t.test(domains.df$self_ind_pct[x==1 & domains.df$rebel==1], domains.df$self_ind_pct[x==1 & domains.df$rebel==0])$estimate[2]})
    TableA7a[i,4]<-lapply(domains.df[domains_top_sig[i]], function(x) {t.test(domains.df$self_ind_pct[x==1 & domains.df$rebel==1], domains.df$self_ind_pct[x==1 & domains.df$rebel==0])$p.value})
  }, error=function(e){})
  TableA7a[i,5]<-lapply(domains.df[domains_top_sig[i]], function(x) {sum(x)})

}

names(TableA7a)<-c("Variable", "Mean Rebel", "Mean Govt", "P-value", "N")
TableA7a$Diff<-TableA7a$`Mean Rebel`- TableA7a$`Mean Govt`
TableA7a<- TableA7a %>% dplyr::select(Variable,  `Mean Rebel`, `Mean Govt`, Diff, `P-value`, N) %>% mutate_if(is.numeric, round, 2)

print(xtable(TableA7a), include.rownames=F, file="TableA7a.tex")


TableA7b<-data.frame(domain_imbal)

for (i in 1:length(domains_top_sig)){
  tryCatch({
  TableA7b[i,2]<-lapply(domains.df[domains_top_sig[i]], function(x) {t.test(domains.df$outgroup_pct[x==1 & domains.df$rebel==1], domains.df$outgroup_pct[x==1 & domains.df$rebel==0])$estimate[1]})
  TableA7b[i,3]<-lapply(domains.df[domains_top_sig[i]], function(x) {t.test(domains.df$outgroup_pct[x==1 & domains.df$rebel==1], domains.df$outgroup_pct[x==1 & domains.df$rebel==0])$estimate[2]})
  TableA7b[i,4]<-lapply(domains.df[domains_top_sig[i]], function(x) {t.test(domains.df$outgroup_pct[x==1 & domains.df$rebel==1], domains.df$outgroup_pct[x==1 & domains.df$rebel==0])$p.value})
    }, error=function(e){})
  TableA7b[i,5]<-lapply(domains.df[domains_top_sig[i]], function(x) {sum(x)})

}

names(TableA7b)<-c("Variable", "Mean Rebel", "Mean Govt", "P-value", "N")
TableA7b$Diff<-TableA7b$`Mean Rebel`- TableA7b$`Mean Govt`
TableA7b<- TableA7b %>% dplyr::select(Variable,  `Mean Rebel`, `Mean Govt`, Diff, `P-value`, N) %>% mutate_if(is.numeric, round, 2)

print(xtable(TableA7b), include.rownames=F, file="TableA7b.tex")

# Table A8: Rate of discrimination against self -------------------------------------------------------------

# select the round 3 variables
r3<-c("R3_disc_self", "R3_disc_self_p", "R3_disc_self_d","r3_always_disc", "r3_once_disc")

r3_2<-c("Ethnicity game", "Partisanship game", "Department game", "All three games", "At least one game")

Table_A8<-data.frame(r3_2)

for (i in 1:length(r3)){
  Table_A8[i,2]<-lapply(labs[r3[i]], function(x) {t.test(x[labs$rebel==1], x[labs$rebel==0])$estimate[1]})
  Table_A8[i,3]<-lapply(labs[r3[i]], function(x) {t.test(x[labs$rebel==1], x[labs$rebel==0])$estimate[2]})
  Table_A8[i,4]<-lapply(labs[r3[i]], function(x) {t.test(x[labs$rebel==1], x[labs$rebel==0])$p.value})
  Table_A8[i,5]<-lapply(labs[r3[i]], function(x) {sum(!is.na(x))})

}


names(Table_A8)<-c("Variable", "Mean Rebel", "Mean Govt", "P-value", "N")
Table_A8$Diff<-Table_A8$`Mean Rebel`- Table_A8$`Mean Govt`
Table_A8<- Table_A8 %>% dplyr::select(Variable,  `Mean Rebel`, `Mean Govt`, Diff, `P-value`, N) %>% mutate_if(is.numeric, round, 2)

print(xtable(Table_A8), include.rownames=F, file="Table_A8.tex")

# Figure A2: Share of discrimination against outgroup once -------------------------------------------------


out_plot_disc<-summarySE(labs, measurevar="discriminate_all", groupvars = c("rebel"), na.rm=T)
out_plot_disc$rebel[out_plot_disc$rebel==1]<-"rebel"
out_plot_disc$rebel[out_plot_disc$rebel==0]<-"government"


ggplot(out_plot_disc, aes(x=rebel, y=100*discriminate_all, fill=rebel)) +
  geom_bar(stat="identity", width=.5, color="black")+
  geom_errorbar(aes(ymin=100*(discriminate_all-ci), ymax=100*(discriminate_all+ci)),
                width = 0.1, position = position_dodge(width = 1/2)) +
  theme_minimal() +
  labs(y = "Discriminate against outgroup \nat least once", x="")+
  scale_x_discrete(labels=c("Leader from \nformer-Govt", "Leader from\nformer-Rebel"))+
  geom_text(aes(label=round(100*out_plot_disc$discriminate_all,2)),
            vjust=-0.5, hjust = -0.1, nudge_x = 0.05, size=7)+
  geom_segment(aes(x=1, y=95, xend=1.1, yend=97), colour="darkgray")+
  geom_segment(aes(x=1.1, y=97, xend=1.9, yend=97), colour="darkgray")+
  geom_segment(aes(x=1.9, y=97, xend=2, yend=95), colour="darkgray")+
  annotate("text", x=1.5, y=100, label="20%p (p<0.01)", size=7)+
  #ylim(0,80)+
  #scale_fill_manual(values=c("#56B4E9", "#0072B2"))+
  scale_fill_manual(values=c("white", "darkgray"))+
  theme(axis.text.x  = element_text(size=20),
        axis.text.y=element_text(size=20),
        axis.title.y = element_text(size=20))

ggsave(filename="FigureA2.pdf", device = cairo_pdf, width = 15, height=8.5, units="in", dpi = 600)


# Table A9: Full regression controlling for individual and org cov --------

stargazer(base_self.lm2, self_ind.lm2, org_ctrls.self, full_reg_self,
          base_out.lm2, out.lm2, org_ctrls_out.out,full_reg_out, style = "apsr",
          type="text",
          column.labels   = c("Percent allocated to self", "Percent allocated to outgroups"),
          column.separate = c(4, 4),
          add.lines = list(c("Observations", paste(nobs(base_self.lm)),paste(nobs(self_ind.lm)),
                             paste(nobs(org_ctrls)),paste(nobs(full_reg2)),
                             paste(nobs(base_out.lm)),paste(nobs(out.lm)),
                             paste(nobs(org_ctrls_out)),paste(nobs(full_reg2))),
                           c("Individual controls", "No", "Yes","No", "Yes","No", "Yes", "No", "Yes"),
                           c("Organization controls", "No", "No","Yes", "Yes","No", "No", "Yes", "Yes")),
covariate.labels = c("Civil society leader from rebel depts", "Sex", "Secondary education or higher",
                     "Northern ethnic group", "Organization ever received international aid",
                     " National headquarters", "Total number of domains",
                     "Created during the war", "Created after the war",
                     "Domain: Humanitarian"," Org mission: information provision"),
out = "TableA9.tex")

# Table A10: regressions with raw data -----------------------------------------------

#remove old regressions
rm(base_self.lm2, self_ind.lm2, org_ctrls.self, full_reg_self,
   base_out.lm2, out.lm2, org_ctrls_out.out,full_reg_out)

base_self.lm<-lm(self_ind~rebel, data=labs)
base_self.lm2<-coeftest(base_self.lm, vcov =  sandwich::vcovCL(base_self.lm, labs$DEPT[is.na(labs$rebel) == FALSE]), df = 20)

self_ind.lm<-lm(self_ind~rebel+SEX+higher_edu+northern, data=labs)
self_ind.lm2<-coeftest(self_ind.lm, vcov =  sandwich::vcovCL(self_ind.lm, labs$DEPT[is.na(labs$rebel) == FALSE  &
                                                                                      is.na(labs$SEX) == FALSE &
                                                                                      is.na(labs$self_ind)==F &
                                                                                      is.na(labs$higher_edu)==F &
                                                                                      is.na(labs$northern)==F]), df = 20)
org_ctrls<-lm(self_ind~rebel+int_fund+natl+total_domains+created_war+created_postwar+DOMAIN_HUMA+ORG_MISSION_INFO_2, data=labs)

org_ctrls.self<-coeftest(org_ctrls, vcov =  sandwich::
                           vcovCL(org_ctrls, labs$DEPT.x[is.na(labs$rebel) == FALSE  &
                                                           is.na(labs$int_fund) == FALSE &
                                                           is.na(labs$self_ind)==F &
                                                           is.na(labs$natl)==F &
                                                           is.na(labs$total_domains)==F &
                                                           is.na(labs$created_postwar)==F &
                                                           is.na(labs$DOMAIN_HUMA)==F &

                                                           is.na(labs$ORG_MISSION_INFO_2)==F &
                                                           is.na(labs$created_war)==F]), df = 20)


base_out.lm<-lm(outgroup~rebel, data=labs)
base_out.lm2<-coeftest(base_out.lm, vcov =  sandwich::vcovCL(base_out.lm, labs$DEPT[is.na(labs$rebel) == FALSE]), df = 20)

out.lm<-lm(outgroup~rebel+SEX+higher_edu+northern, data=labs)
out.lm2<-coeftest(out.lm, vcov =  sandwich::vcovCL(out.lm, labs$DEPT[is.na(labs$rebel) == FALSE  &
                                                                       is.na(labs$SEX) == FALSE &
                                                                       is.na(labs$outgroup)==F &
                                                                       is.na(labs$higher_edu)==F &
                                                                       is.na(labs$northern)==F]), df = 20)
org_ctrls_out<-lm(outgroup~rebel+int_fund+natl+total_domains+created_war+created_postwar+DOMAIN_HUMA+ORG_MISSION_INFO_2, data=labs)
org_ctrls_out.out<-coeftest(org_ctrls_out, vcov =  sandwich::
                              vcovCL(org_ctrls_out, labs$DEPT.x[is.na(labs$rebel) == FALSE  &
                                                                  is.na(labs$int_fund) == FALSE &
                                                                  is.na(labs$outgroup)==F &
                                                                  is.na(labs$natl)==F &
                                                                  is.na(labs$total_domains)==F &
                                                                  is.na(labs$created_postwar)==F &
                                                                  is.na(labs$DOMAIN_HUMA)==F &

                                                                  is.na(labs$ORG_MISSION_INFO_2)==F &
                                                                  is.na(labs$created_war)==F]), df = 20)



full_reg<-lm(self_ind~rebel+int_fund+natl+total_domains+created_war+created_postwar+DOMAIN_HUMA+ORG_MISSION_INFO_2+SEX+higher_edu+northern, data=labs)
full_reg_self<-coeftest(full_reg, vcov =  sandwich::
                          vcovCL(full_reg, labs$DEPT.x[is.na(labs$rebel) == FALSE  &
                                                         is.na(labs$int_fund) == FALSE &
                                                         is.na(labs$self_ind)==F &
                                                         is.na(labs$natl)==F &
                                                         is.na(labs$total_domains)==F &
                                                         is.na(labs$created_postwar)==F &
                                                         is.na(labs$DOMAIN_HUMA)==F &
                                                         is.na(labs$SEX) == FALSE &
                                                         is.na(labs$higher_edu)==F &
                                                         is.na(labs$northern)==F &
                                                         is.na(labs$ORG_MISSION_INFO_2)==F &
                                                         is.na(labs$created_war)==F]), df = 20)

full_reg2<-lm(outgroup~rebel+int_fund+natl+total_domains+created_war+created_postwar+DOMAIN_HUMA+ORG_MISSION_INFO_2+SEX+higher_edu+northern, data=labs)
full_reg_out<-coeftest(full_reg2, vcov =  sandwich::
                         vcovCL(full_reg2, labs$DEPT.x[is.na(labs$rebel) == FALSE  &
                                                         is.na(labs$int_fund) == FALSE &
                                                         is.na(labs$outgroup)==F &
                                                         is.na(labs$natl)==F &
                                                         is.na(labs$total_domains)==F &
                                                         is.na(labs$created_postwar)==F &
                                                         is.na(labs$DOMAIN_HUMA)==F &
                                                         is.na(labs$SEX) == FALSE &
                                                         is.na(labs$higher_edu)==F &
                                                         is.na(labs$northern)==F &
                                                         is.na(labs$ORG_MISSION_INFO_2)==F &
                                                         is.na(labs$created_war)==F]), df = 20)





require(stargazer)
stargazer(base_self.lm2, self_ind.lm2, org_ctrls.self, full_reg_self,
          base_out.lm2, out.lm2, org_ctrls_out.out,full_reg_out, style = "apsr",
          type="text",
          column.labels   = c("Percent allocated to self", "Percent allocated to outgroups"),
          column.separate = c(4, 4),
          add.lines = list(c("Observations", paste(nobs(base_self.lm)),paste(nobs(self_ind.lm)),
                             paste(nobs(org_ctrls)),paste(nobs(full_reg2)),
                             paste(nobs(base_out.lm)),paste(nobs(out.lm)),
                             paste(nobs(org_ctrls_out)),paste(nobs(full_reg2))),
                           c("Individual controls", "No", "Yes","No", "Yes","No", "Yes", "No", "Yes"),
                           c("Organization controls", "No", "No","Yes", "Yes","No", "No", "Yes", "Yes")),
          keep = c("rebel"),
          covariate.labels = "Civil society leader from rebel depts",
          out = "TableA10.tex")
# Table A11: Outgroup analysis as share of amount that could be allocated to the outgroup ---------------------------------------------------


#remove old regressions
rm(base_self.lm2, self_ind.lm2, org_ctrls.self, full_reg_self,
   base_out.lm2, out.lm2, org_ctrls_out.out,full_reg_out)

base_out.lm<-lm(outgroup_pct_total~rebel, data=labs)
base_out.lm2<-coeftest(base_out.lm, vcov =  sandwich::vcovCL(base_out.lm, labs$DEPT[is.na(labs$rebel) == FALSE]), df = 20)

out.lm<-lm(outgroup_pct_total~rebel+SEX+higher_edu+northern, data=labs)
out.lm2<-coeftest(out.lm, vcov =  sandwich::vcovCL(out.lm, labs$DEPT[is.na(labs$rebel) == FALSE  &
                                                                       is.na(labs$SEX) == FALSE &
                                                                       is.na(labs$outgroup_pct_total)==F &
                                                                       is.na(labs$higher_edu)==F &
                                                                       is.na(labs$northern)==F]), df = 20)
org_ctrls_out<-lm(outgroup_pct_total~rebel+int_fund+natl+total_domains+created_war+created_postwar+DOMAIN_HUMA+ORG_MISSION_INFO_2, data=labs)
org_ctrls_out.out<-coeftest(org_ctrls_out, vcov =  sandwich::
                              vcovCL(org_ctrls_out, labs$DEPT.x[is.na(labs$rebel) == FALSE  &
                                                                  is.na(labs$int_fund) == FALSE &
                                                                  is.na(labs$outgroup_pct_total)==F &
                                                                  is.na(labs$natl)==F &
                                                                  is.na(labs$total_domains)==F &
                                                                  is.na(labs$created_postwar)==F &
                                                                  is.na(labs$DOMAIN_HUMA)==F &

                                                                  is.na(labs$ORG_MISSION_INFO_2)==F &
                                                                  is.na(labs$created_war)==F]), df = 20)



full_reg2<-lm(outgroup_pct_total~rebel+int_fund+natl+total_domains+created_war+created_postwar+DOMAIN_HUMA+ORG_MISSION_INFO_2+SEX+higher_edu+northern, data=labs)
full_reg_out<-coeftest(full_reg2, vcov =  sandwich::
                         vcovCL(full_reg2, labs$DEPT.x[is.na(labs$rebel) == FALSE  &
                                                         is.na(labs$int_fund) == FALSE &
                                                         is.na(labs$outgroup_pct_total)==F &
                                                         is.na(labs$natl)==F &
                                                         is.na(labs$total_domains)==F &
                                                         is.na(labs$created_postwar)==F &
                                                         is.na(labs$DOMAIN_HUMA)==F &
                                                         is.na(labs$SEX) == FALSE &
                                                         is.na(labs$higher_edu)==F &
                                                         is.na(labs$northern)==F &
                                                         is.na(labs$ORG_MISSION_INFO_2)==F &
                                                         is.na(labs$created_war)==F]), df = 20)



require(stargazer)
stargazer(base_out.lm2, out.lm2, org_ctrls_out.out,full_reg_out, style = "apsr",
          type="text",
          column.labels   = c("Percent allocated to outgroup, out of remaining total"),
          column.separate = c(4),
          add.lines = list(c("Observations", paste(nobs(base_self.lm)),paste(nobs(self_ind.lm)),
                             paste(nobs(org_ctrls)),paste(nobs(full_reg2)),
                             paste(nobs(base_out.lm)),paste(nobs(out.lm)),
                             paste(nobs(org_ctrls_out)),paste(nobs(full_reg2))),
                           c("Individual controls", "No", "Yes","No", "Yes","No", "Yes", "No", "Yes"),
                           c("Organization controls", "No", "No","Yes", "Yes","No", "No", "Yes", "Yes")),
          keep = c("rebel"),
          covariate.labels = "Civil society leader from rebel depts",
          out = "TableA11.tex")



# Table A12: Cacao production ----------------------------------------------

# to examine cacao production, we exclude the outlying departments and examine cacao production to make sure there is not imbalance on cacao production
zdc<- nielsen_match_data %>% filter(dept!="Mankono") %>%
  filter(dept!="Niakaramandougou") %>%
  filter(dept!="Sassandra") %>%
  filter(dept!="San Pedro") %>%
  filter(dept!="Tabou")

Table_A12<-data.frame(variable="Cacao production")

Table_A12[1,2]<-t.test(zdc$cacao_production[zdc$rebel==1], zdc$cacao_production[zdc$rebel==0])$estimate[1]
Table_A12[1,3]<-t.test(zdc$cacao_production[zdc$rebel==1], zdc$cacao_production[zdc$rebel==0])$estimate[2]
Table_A12[1,4]<-t.test(zdc$cacao_production[zdc$rebel==1], zdc$cacao_production[zdc$rebel==0])$p.value
Table_A12[1,5]<-nrow(zdc)


names(Table_A12)<-c("Variable", "Mean Rebel", "Mean Govt", "P-value", "N")
Table_A12$Diff<-Table_A12$`Mean Rebel`- Table_A12$`Mean Govt`
Table_A12<- Table_A12 %>% dplyr::select(Variable,  `Mean Rebel`, `Mean Govt`, Diff, `P-value`, N) %>% mutate_if(is.numeric, round, 2)

print(xtable(Table_A12), include.rownames=F, file="Table_A12.tex")



# Table A13: Analysis excluding outlying depts ----------------------------

labs_zdc<- labs %>% filter(zdc==1)

#remove old regressions
rm(base_self.lm2, self_ind.lm2, org_ctrls.self, full_reg_self,
   base_out.lm2, out.lm2, org_ctrls_out.out,full_reg_out)

base_self.lm<-lm(self_ind_pct~rebel, data=labs_zdc)
base_self.lm2<-coeftest(base_self.lm, vcov =  sandwich::vcovCL(base_self.lm, labs_zdc$DEPT[is.na(labs_zdc$rebel) == FALSE]), df = 15)

self_ind.lm<-lm(self_ind_pct~rebel+SEX+higher_edu+northern, data=labs_zdc)
self_ind.lm2<-coeftest(self_ind.lm, vcov =  sandwich::vcovCL(self_ind.lm, labs_zdc$DEPT[is.na(labs_zdc$rebel) == FALSE  &
                                                                                          is.na(labs_zdc$SEX) == FALSE &
                                                                                          is.na(labs_zdc$self_ind_pct)==F &
                                                                                          is.na(labs_zdc$higher_edu)==F &
                                                                                          is.na(labs_zdc$northern)==F]), df = 15)
org_ctrls<-lm(self_ind_pct~rebel+int_fund+natl+total_domains+created_war+created_postwar+DOMAIN_HUMA+ORG_MISSION_INFO_2, data=labs_zdc)

org_ctrls.self<-coeftest(org_ctrls, vcov =  sandwich::
                           vcovCL(org_ctrls, labs_zdc$DEPT.x[is.na(labs_zdc$rebel) == FALSE  &
                                                               is.na(labs_zdc$int_fund) == FALSE &
                                                               is.na(labs_zdc$self_ind_pct)==F &
                                                               is.na(labs_zdc$natl)==F &
                                                               is.na(labs_zdc$total_domains)==F &
                                                               is.na(labs_zdc$created_postwar)==F &
                                                               is.na(labs_zdc$DOMAIN_HUMA)==F &
                                                               is.na(labs_zdc$ORG_MISSION_INFO_2)==F &
                                                               is.na(labs_zdc$created_war)==F]), df = 15)


base_out.lm<-lm(outgroup_pct~rebel, data=labs_zdc)
base_out.lm2<-coeftest(base_out.lm, vcov =  sandwich::vcovCL(base_out.lm, labs_zdc$DEPT[is.na(labs_zdc$rebel) == FALSE]), df = 15)

out.lm<-lm(outgroup_pct~rebel+SEX+higher_edu+northern, data=labs_zdc)
out.lm2<-coeftest(out.lm, vcov =  sandwich::vcovCL(out.lm, labs_zdc$DEPT[is.na(labs_zdc$rebel) == FALSE  &
                                                                           is.na(labs_zdc$SEX) == FALSE &
                                                                           is.na(labs_zdc$outgroup_pct)==F &
                                                                           is.na(labs_zdc$higher_edu)==F &
                                                                           is.na(labs_zdc$northern)==F]), df = 15)
org_ctrls_out<-lm(outgroup_pct~rebel+int_fund+natl+total_domains+created_war+created_postwar+DOMAIN_HUMA+ORG_MISSION_INFO_2, data=labs_zdc)
org_ctrls_out.out<-coeftest(org_ctrls_out, vcov =  sandwich::
                              vcovCL(org_ctrls_out, labs_zdc$DEPT.x[is.na(labs_zdc$rebel) == FALSE  &
                                                                      is.na(labs_zdc$int_fund) == FALSE &
                                                                      is.na(labs_zdc$outgroup_pct)==F &
                                                                      is.na(labs_zdc$natl)==F &
                                                                      is.na(labs_zdc$total_domains)==F &
                                                                      is.na(labs_zdc$created_postwar)==F &
                                                                      is.na(labs_zdc$DOMAIN_HUMA)==F &

                                                                      is.na(labs_zdc$ORG_MISSION_INFO_2)==F &
                                                                      is.na(labs_zdc$created_war)==F]), df = 15)



full_reg<-lm(self_ind_pct~rebel+int_fund+natl+total_domains+created_war+created_postwar+DOMAIN_HUMA+ORG_MISSION_INFO_2+SEX+higher_edu+northern, data=labs_zdc)
full_reg_self<-coeftest(full_reg, vcov =  sandwich::
                          vcovCL(full_reg, labs_zdc$DEPT.x[is.na(labs_zdc$rebel) == FALSE  &
                                                             is.na(labs_zdc$int_fund) == FALSE &
                                                             is.na(labs_zdc$self_ind_pct)==F &
                                                             is.na(labs_zdc$natl)==F &
                                                             is.na(labs_zdc$total_domains)==F &
                                                             is.na(labs_zdc$created_postwar)==F &
                                                             is.na(labs_zdc$DOMAIN_HUMA)==F &
                                                             is.na(labs_zdc$SEX) == FALSE &
                                                             is.na(labs_zdc$higher_edu)==F &
                                                             is.na(labs_zdc$northern)==F &
                                                             is.na(labs_zdc$ORG_MISSION_INFO_2)==F &
                                                             is.na(labs_zdc$created_war)==F]), df = 15)

full_reg2<-lm(outgroup_pct~rebel+int_fund+natl+total_domains+created_war+created_postwar+DOMAIN_HUMA+ORG_MISSION_INFO_2+SEX+higher_edu+northern, data=labs_zdc)
full_reg_out<-coeftest(full_reg2, vcov =  sandwich::
                         vcovCL(full_reg2, labs_zdc$DEPT.x[is.na(labs_zdc$rebel) == FALSE  &
                                                             is.na(labs_zdc$int_fund) == FALSE &
                                                             is.na(labs_zdc$outgroup_pct)==F &
                                                             is.na(labs_zdc$natl)==F &
                                                             is.na(labs_zdc$total_domains)==F &
                                                             is.na(labs_zdc$created_postwar)==F &
                                                             is.na(labs_zdc$DOMAIN_HUMA)==F &
                                                             is.na(labs_zdc$SEX) == FALSE &
                                                             is.na(labs_zdc$higher_edu)==F &
                                                             is.na(labs_zdc$northern)==F &
                                                             is.na(labs_zdc$ORG_MISSION_INFO_2)==F &
                                                             is.na(labs_zdc$created_war)==F]), df = 15)





require(stargazer)
stargazer(base_self.lm2, self_ind.lm2, org_ctrls.self, full_reg_self,
          base_out.lm2, out.lm2, org_ctrls_out.out,full_reg_out, style = "apsr",
          type="text",
          column.labels   = c("Percent allocated to self", "Percent allocated to outgroups"),
          column.separate = c(4, 4),
          add.lines = list(c("Observations", paste(nobs(base_self.lm)),paste(nobs(self_ind.lm)),
                             paste(nobs(org_ctrls)),paste(nobs(full_reg2)),
                             paste(nobs(base_out.lm)),paste(nobs(out.lm)),
                             paste(nobs(org_ctrls_out)),paste(nobs(full_reg2))),
                           c("Individual controls", "No", "Yes","No", "Yes","No", "Yes", "No", "Yes"),
                           c("Organization controls", "No", "No","Yes", "Yes","No", "No", "Yes", "Yes")),
          keep = c("rebel"),
          covariate.labels = "Civil society leader from rebel depts",
          out = "TableA13.tex")


# Figure A3: Cacao production ---------------------------------------------


dept_centroids_cases_cacao<-dept_centroids %>% filter(dept_name %in% nielsen_cases)

ggplot(civ.df)+
  aes(x=long, y=lat) +
  geom_polygon(aes(group=group, fill=factor(civ.df$cacao)), color="gray")+
  coord_map("mercator")+
  scale_fill_manual(values=c("lightgreen", "limegreen", "forestgreen","White"),
                    labels=c("Low production", "Average production", "High production"),
                    name="Côte d'Ivoire, cacao production", na.translate = F)+
  geom_text(data=dept_centroids_cases_cacao, aes(label = dept_name, x = long, y = lat), color="black", size=3)+
  geom_path(data=case_map, aes(x=long, y=lat, group=group), color="black", size=1)+
  geom_line(data=zdc_lower, col="black",linetype = 2, size=1)+
  geom_line(data=zdc_upper, col="black",linetype = 2, size=1)+
  theme(legend.text = element_text(size=20), legend.title = element_text(size=20))+
  theme_classic()

ggsave(filename="FigA3.pdf", device = cairo_pdf, width = 15, height=8.5, units="in", dpi = 600)


# Table A14: Controlling for private funding ------------------------------

#remove old regressions
rm(base_self.lm2, self_ind.lm2, org_ctrls.self, full_reg_self,
   base_out.lm2, out.lm2, org_ctrls_out.out,full_reg_out)

base_self.lm_c<-lm(self_ind_pct~rebel+fund_private, data=labs)
base_self.lm_c2<-coeftest(base_self.lm_c, vcov =  sandwich::vcovCL(base_self.lm_c, labs$DEPT[is.na(labs$rebel) == FALSE &
                                                                                             is.na(labs$fund_private) == FALSE]), df = 20)

base_out.lm_c<-lm(outgroup_pct~rebel+fund_private, data=labs)
base_out.lm_c2<-coeftest(base_out.lm_c, vcov =  sandwich::vcovCL(base_out.lm_c, labs$DEPT[is.na(labs$rebel) == FALSE &
                                                                                            is.na(labs$outgroup_pct) == FALSE &
                                                                                            is.na(labs$fund_private) == FALSE]), df = 20)

stargazer(base_self.lm_c2, base_out.lm_c2, style = "apsr",
          type="text",
          column.labels   = c("Percent allocated to self", "Percent allocated to outgroups"),
          add.lines = list(c("Observations", paste(nobs(base_self.lm_c2)),paste(nobs(base_out.lm_c2)))),
          keep = c("rebel", "fund_private"),
          covariate.labels = "Civil society leader from rebel depts",
          out = "TableA14.tex")



# Table A15: Violence Analysis ----------------------------------------

#remove old regressions
rm(base_self.lm2, self_ind.lm2, org_ctrls.self, full_reg_self,
   base_out.lm2, out.lm2, org_ctrls_out.out,full_reg_out)

base_self.lm<-lm(self_ind_pct~rebel+violence2, data=labs)
base_self.lm2<-coeftest(base_self.lm, vcov =  sandwich::vcovCL(base_self.lm, labs$DEPT[is.na(labs$rebel) == FALSE]), df = 20)

self_ind.lm<-lm(self_ind_pct~rebel+violence2+SEX+higher_edu+northern, data=labs)
self_ind.lm2<-coeftest(self_ind.lm, vcov =  sandwich::vcovCL(self_ind.lm, labs$DEPT[is.na(labs$rebel) == FALSE  &
                                                                                      is.na(labs$SEX) == FALSE &
                                                                                      is.na(labs$self_ind_pct)==F &
                                                                                      is.na(labs$higher_edu)==F &
                                                                                      is.na(labs$violence2)==F &
                                                                                      is.na(labs$northern)==F]), df = 20)
org_ctrls<-lm(self_ind_pct~rebel+int_fund+natl+total_domains+created_war+created_postwar+DOMAIN_HUMA+ORG_MISSION_INFO_2, data=labs)

org_ctrls.self<-coeftest(org_ctrls, vcov =  sandwich::
                           vcovCL(org_ctrls, labs$DEPT.x[is.na(labs$rebel) == FALSE  &
                                                           is.na(labs$int_fund) == FALSE &
                                                           is.na(labs$self_ind_pct)==F &
                                                           is.na(labs$natl)==F &
                                                           is.na(labs$total_domains)==F &
                                                           is.na(labs$created_postwar)==F &
                                                           is.na(labs$DOMAIN_HUMA)==F &

                                                           is.na(labs$ORG_MISSION_INFO_2)==F &
                                                           is.na(labs$created_war)==F]), df = 20)


base_out.lm<-lm(outgroup_pct~rebel+violence2, data=labs)
base_out.lm2<-coeftest(base_out.lm, vcov =  sandwich::vcovCL(base_out.lm, labs$DEPT[is.na(labs$rebel) == FALSE]), df = 20)

out.lm<-lm(outgroup_pct~rebel+violence2+SEX+higher_edu+northern, data=labs)
out.lm2<-coeftest(out.lm, vcov =  sandwich::vcovCL(out.lm, labs$DEPT[is.na(labs$rebel) == FALSE  &
                                                                       is.na(labs$SEX) == FALSE &
                                                                       is.na(labs$outgroup_pct)==F &
                                                                       is.na(labs$violence2)==F &
                                                                       is.na(labs$higher_edu)==F &
                                                                       is.na(labs$northern)==F]), df = 20)
org_ctrls_out<-lm(outgroup_pct~rebel+violence2+int_fund+natl+total_domains+created_war+created_postwar+DOMAIN_HUMA+ORG_MISSION_INFO_2, data=labs)
org_ctrls_out.out<-coeftest(org_ctrls_out, vcov =  sandwich::
                              vcovCL(org_ctrls_out, labs$DEPT.x[is.na(labs$rebel) == FALSE  &
                                                                  is.na(labs$int_fund) == FALSE &
                                                                  is.na(labs$outgroup_pct)==F &
                                                                  is.na(labs$natl)==F &
                                                                  is.na(labs$total_domains)==F &
                                                                  is.na(labs$violence2)==F &
                                                                  is.na(labs$created_postwar)==F &
                                                                  is.na(labs$DOMAIN_HUMA)==F &

                                                                  is.na(labs$ORG_MISSION_INFO_2)==F &
                                                                  is.na(labs$created_war)==F]), df = 20)



full_reg<-lm(self_ind_pct~rebel+violence2+int_fund+natl+total_domains+created_war+created_postwar+DOMAIN_HUMA+ORG_MISSION_INFO_2+SEX+higher_edu+northern, data=labs)
full_reg_self<-coeftest(full_reg, vcov =  sandwich::
                          vcovCL(full_reg, labs$DEPT.x[is.na(labs$rebel) == FALSE  &
                                                         is.na(labs$int_fund) == FALSE &
                                                         is.na(labs$self_ind_pct)==F &
                                                         is.na(labs$natl)==F &
                                                         is.na(labs$total_domains)==F &
                                                         is.na(labs$created_postwar)==F &
                                                         is.na(labs$violence2)==F &
                                                         is.na(labs$DOMAIN_HUMA)==F &
                                                         is.na(labs$SEX) == FALSE &
                                                         is.na(labs$higher_edu)==F &
                                                         is.na(labs$northern)==F &
                                                         is.na(labs$ORG_MISSION_INFO_2)==F &
                                                         is.na(labs$created_war)==F]), df = 20)

full_reg2<-lm(outgroup_pct~rebel+violence2+int_fund+natl+total_domains+created_war+created_postwar+DOMAIN_HUMA+ORG_MISSION_INFO_2+SEX+higher_edu+northern, data=labs)
full_reg_out<-coeftest(full_reg2, vcov =  sandwich::
                         vcovCL(full_reg2, labs$DEPT.x[is.na(labs$rebel) == FALSE  &
                                                         is.na(labs$int_fund) == FALSE &
                                                         is.na(labs$outgroup_pct)==F &
                                                         is.na(labs$natl)==F &
                                                         is.na(labs$total_domains)==F &
                                                         is.na(labs$created_postwar)==F &
                                                         is.na(labs$violence2)==F &
                                                         is.na(labs$DOMAIN_HUMA)==F &
                                                         is.na(labs$SEX) == FALSE &
                                                         is.na(labs$higher_edu)==F &
                                                         is.na(labs$northern)==F &
                                                         is.na(labs$ORG_MISSION_INFO_2)==F &
                                                         is.na(labs$created_war)==F]), df = 20)





require(stargazer)
stargazer(base_self.lm2, base_out.lm2, style = "apsr",
          type="text",
          column.labels   = c("Percent allocated to self", "Percent allocated to outgroups"),
          column.separate = c(1, 1),
          add.lines = list(c("Observations", paste(nobs(base_self.lm2)),
                             paste(nobs(base_out.lm2)))),
          keep = c("rebel", "violence2"),
          covariate.labels = c("Civil society leader from rebel depts", "violence exposure"),
          out = "TableA15.tex")


# Table A16: Analysis including ministry partnerships ------------------------------------------------------

# make sure not to remove old regressions
rm(base_self.lm2, self_ind.lm2, org_ctrls.self, full_reg_self,
   base_out.lm2, out.lm2, org_ctrls_out.out,full_reg_out)

base_self.lm<-lm(self_ind_pct~rebel+min_part, data=labs)
base_self.lm2<-coeftest(base_self.lm, vcov =  sandwich::vcovCL(base_self.lm, labs$DEPT[is.na(labs$rebel) == FALSE]), df = 20)

self_ind.lm<-lm(self_ind_pct~rebel+min_part+SEX+higher_edu+northern, data=labs)
self_ind.lm2<-coeftest(self_ind.lm, vcov =  sandwich::vcovCL(self_ind.lm, labs$DEPT[is.na(labs$rebel) == FALSE  &
                                                                                      is.na(labs$SEX) == FALSE &
                                                                                      is.na(labs$self_ind_pct)==F &
                                                                                      is.na(labs$higher_edu)==F &
                                                                                      is.na(labs$min_part)==F &
                                                                                      is.na(labs$northern)==F]), df = 20)
org_ctrls<-lm(self_ind_pct~rebel+int_fund+natl+total_domains+created_war+created_postwar+DOMAIN_HUMA+ORG_MISSION_INFO_2, data=labs)

org_ctrls.self<-coeftest(org_ctrls, vcov =  sandwich::
                           vcovCL(org_ctrls, labs$DEPT.x[is.na(labs$rebel) == FALSE  &
                                                           is.na(labs$int_fund) == FALSE &
                                                           is.na(labs$self_ind_pct)==F &
                                                           is.na(labs$natl)==F &
                                                           is.na(labs$total_domains)==F &
                                                           is.na(labs$created_postwar)==F &
                                                           is.na(labs$DOMAIN_HUMA)==F &

                                                           is.na(labs$ORG_MISSION_INFO_2)==F &
                                                           is.na(labs$created_war)==F]), df = 20)


base_out.lm<-lm(outgroup_pct~rebel+min_part, data=labs)
base_out.lm2<-coeftest(base_out.lm, vcov =  sandwich::vcovCL(base_out.lm, labs$DEPT[is.na(labs$rebel) == FALSE]), df = 20)

out.lm<-lm(outgroup_pct~rebel+min_part+SEX+higher_edu+northern, data=labs)
out.lm2<-coeftest(out.lm, vcov =  sandwich::vcovCL(out.lm, labs$DEPT[is.na(labs$rebel) == FALSE  &
                                                                       is.na(labs$SEX) == FALSE &
                                                                       is.na(labs$outgroup_pct)==F &
                                                                       is.na(labs$min_part)==F &
                                                                       is.na(labs$higher_edu)==F &
                                                                       is.na(labs$northern)==F]), df = 20)
org_ctrls_out<-lm(outgroup_pct~rebel+min_part+int_fund+natl+total_domains+created_war+created_postwar+DOMAIN_HUMA+ORG_MISSION_INFO_2, data=labs)
org_ctrls_out.out<-coeftest(org_ctrls_out, vcov =  sandwich::
                              vcovCL(org_ctrls_out, labs$DEPT.x[is.na(labs$rebel) == FALSE  &
                                                                  is.na(labs$int_fund) == FALSE &
                                                                  is.na(labs$outgroup_pct)==F &
                                                                  is.na(labs$natl)==F &
                                                                  is.na(labs$total_domains)==F &
                                                                  is.na(labs$min_part)==F &
                                                                  is.na(labs$created_postwar)==F &
                                                                  is.na(labs$DOMAIN_HUMA)==F &

                                                                  is.na(labs$ORG_MISSION_INFO_2)==F &
                                                                  is.na(labs$created_war)==F]), df = 20)



full_reg<-lm(self_ind_pct~rebel+min_part+int_fund+natl+total_domains+created_war+created_postwar+DOMAIN_HUMA+ORG_MISSION_INFO_2+SEX+higher_edu+northern, data=labs)
full_reg_self<-coeftest(full_reg, vcov =  sandwich::
                          vcovCL(full_reg, labs$DEPT.x[is.na(labs$rebel) == FALSE  &
                                                         is.na(labs$int_fund) == FALSE &
                                                         is.na(labs$self_ind_pct)==F &
                                                         is.na(labs$natl)==F &
                                                         is.na(labs$total_domains)==F &
                                                         is.na(labs$created_postwar)==F &
                                                         is.na(labs$min_part)==F &
                                                         is.na(labs$DOMAIN_HUMA)==F &
                                                         is.na(labs$SEX) == FALSE &
                                                         is.na(labs$higher_edu)==F &
                                                         is.na(labs$northern)==F &
                                                         is.na(labs$ORG_MISSION_INFO_2)==F &
                                                         is.na(labs$created_war)==F]), df = 20)

full_reg2<-lm(outgroup_pct~rebel+min_part+int_fund+natl+total_domains+created_war+created_postwar+DOMAIN_HUMA+ORG_MISSION_INFO_2+SEX+higher_edu+northern, data=labs)
full_reg_out<-coeftest(full_reg2, vcov =  sandwich::
                         vcovCL(full_reg2, labs$DEPT.x[is.na(labs$rebel) == FALSE  &
                                                         is.na(labs$int_fund) == FALSE &
                                                         is.na(labs$outgroup_pct)==F &
                                                         is.na(labs$natl)==F &
                                                         is.na(labs$total_domains)==F &
                                                         is.na(labs$created_postwar)==F &
                                                         is.na(labs$min_part)==F &
                                                         is.na(labs$DOMAIN_HUMA)==F &
                                                         is.na(labs$SEX) == FALSE &
                                                         is.na(labs$higher_edu)==F &
                                                         is.na(labs$northern)==F &
                                                         is.na(labs$ORG_MISSION_INFO_2)==F &
                                                         is.na(labs$created_war)==F]), df = 20)





require(stargazer)
stargazer(base_self.lm2, base_out.lm2, style = "apsr",
          type="text",
          column.labels   = c("Percent allocated to self", "Percent allocated to outgroups"),
          column.separate = c(1, 1),
          add.lines = list(c("Observations", paste(nobs(base_self.lm2)),
                             paste(nobs(base_out.lm2)))),
          keep = c("rebel", "min_part"),
          covariate.labels = c("Civil society leader from rebel depts", "Ministry partnerships"),
          out = "TableA16.tex")


# Tables A17: Wartime leader analysis - self --------------------------------

labs_wl<- labs %>% filter(wartime_leaders==1)
labs_pl<- labs %>% filter(postwar_leader==1)


TableA17<-data.frame(variable=c("wartime leader", "postwar leader"))
TableA17[1,2]<- t.test(labs_wl$self_ind_pct[labs_wl$rebel==1], labs_wl$self_ind_pct[labs_wl$rebel==0])$estimate[1]
TableA17[2,2]<- t.test(labs_pl$self_ind_pct[labs_pl$rebel==1], labs_pl$self_ind_pct[labs_pl$rebel==0])$estimate[1]
TableA17[1,3]<- t.test(labs_wl$self_ind_pct[labs_wl$rebel==1], labs_wl$self_ind_pct[labs_wl$rebel==0])$estimate[2]
TableA17[2,3]<- t.test(labs_pl$self_ind_pct[labs_pl$rebel==1], labs_pl$self_ind_pct[labs_pl$rebel==0])$estimate[2]
TableA17[1,4]<- t.test(labs_wl$self_ind_pct[labs_wl$rebel==1], labs_wl$self_ind_pct[labs_wl$rebel==0])$p.value
TableA17[2,4]<- t.test(labs_pl$self_ind_pct[labs_pl$rebel==1], labs_pl$self_ind_pct[labs_pl$rebel==0])$p.value
TableA17[1,5]<- nrow(labs_wl)
TableA17[2,5]<- nrow(labs_pl)


names(TableA17)<-c("Variable", "Mean Rebel", "Mean Govt", "P-value", "N")
TableA17<- TableA17 %>% mutate_if(is.numeric, round, 2)

print(xtable(TableA17), include.rownames=F, file="TableA17.tex")

# Tables A18: Wartime leader analysis - outgroup --------------------------------

TableA18<-data.frame(variable=c("wartime leader", "postwar leader"))
TableA18[1,2]<- t.test(labs_wl$outgroup_pct[labs_wl$rebel==1], labs_wl$outgroup_pct[labs_wl$rebel==0])$estimate[1]
TableA18[2,2]<- t.test(labs_pl$outgroup_pct[labs_pl$rebel==1], labs_pl$outgroup_pct[labs_pl$rebel==0])$estimate[1]
TableA18[1,3]<- t.test(labs_wl$outgroup_pct[labs_wl$rebel==1], labs_wl$outgroup_pct[labs_wl$rebel==0])$estimate[2]
TableA18[2,3]<- t.test(labs_pl$outgroup_pct[labs_pl$rebel==1], labs_pl$outgroup_pct[labs_pl$rebel==0])$estimate[2]
TableA18[1,4]<- t.test(labs_wl$outgroup_pct[labs_wl$rebel==1], labs_wl$outgroup_pct[labs_wl$rebel==0])$p.value
TableA18[2,4]<- t.test(labs_pl$outgroup_pct[labs_pl$rebel==1], labs_pl$outgroup_pct[labs_pl$rebel==0])$p.value
TableA18[1,5]<- nrow(labs_wl)
TableA18[2,5]<- nrow(labs_pl)

names(TableA18)<-c("Variable", "Mean Rebel", "Mean Govt", "P-value", "N")
TableA18<- TableA18 %>% mutate_if(is.numeric, round, 2)

print(xtable(TableA18), include.rownames=F, file="TableA18.tex")

# Table A19: Wartime leaders only regression----------------------------------------------------

#remove old regressions
rm(base_self.lm2, self_ind.lm2, org_ctrls.self, full_reg_self,
   base_out.lm2, out.lm2, org_ctrls_out.out,full_reg_out)

base_self.lm<-lm(self_ind_pct~rebel, data=labs_wl)
base_self.lm2<-coeftest(base_self.lm, vcov =  sandwich::vcovCL(base_self.lm, labs_wl$DEPT[is.na(labs_wl$rebel) == FALSE]), df = 20)

self_ind.lm<-lm(self_ind_pct~rebel+SEX+higher_edu+northern, data=labs_wl)
self_ind.lm2<-coeftest(self_ind.lm, vcov =  sandwich::vcovCL(self_ind.lm, labs_wl$DEPT[is.na(labs_wl$rebel) == FALSE  &
                                                                                      is.na(labs_wl$SEX) == FALSE &
                                                                                      is.na(labs_wl$self_ind_pct)==F &
                                                                                      is.na(labs_wl$higher_edu)==F &
                                                                                      is.na(labs_wl$northern)==F]), df = 20)
org_ctrls<-lm(self_ind_pct~rebel+int_fund+natl+total_domains+DOMAIN_HUMA+ORG_MISSION_INFO_2, data=labs_wl)

org_ctrls.self<-coeftest(org_ctrls, vcov =  sandwich::
                           vcovCL(org_ctrls, labs_wl$DEPT.x[is.na(labs_wl$rebel) == FALSE  &
                                                           is.na(labs_wl$int_fund) == FALSE &
                                                           is.na(labs_wl$self_ind_pct)==F &
                                                           is.na(labs_wl$natl)==F &
                                                           is.na(labs_wl$total_domains)==F &
                                                           is.na(labs_wl$DOMAIN_HUMA)==F &

                                                           is.na(labs_wl$ORG_MISSION_INFO_2)==F]), df = 20)


base_out.lm<-lm(outgroup_pct~rebel+violence2, data=labs_wl)
base_out.lm2<-coeftest(base_out.lm, vcov =  sandwich::vcovCL(base_out.lm, labs_wl$DEPT[is.na(labs_wl$rebel) == FALSE]), df = 20)

out.lm<-lm(outgroup_pct~rebel+SEX+higher_edu+northern, data=labs_wl)
out.lm2<-coeftest(out.lm, vcov =  sandwich::vcovCL(out.lm, labs_wl$DEPT[is.na(labs_wl$rebel) == FALSE  &
                                                                       is.na(labs_wl$SEX) == FALSE &
                                                                       is.na(labs_wl$outgroup_pct)==F &
                                                                       is.na(labs_wl$higher_edu)==F &
                                                                       is.na(labs_wl$northern)==F]), df = 20)
org_ctrls_out<-lm(outgroup_pct~rebel+int_fund+natl+total_domains+DOMAIN_HUMA+ORG_MISSION_INFO_2, data=labs_wl)
org_ctrls_out.out<-coeftest(org_ctrls_out, vcov =  sandwich::
                              vcovCL(org_ctrls_out, labs_wl$DEPT.x[is.na(labs_wl$rebel) == FALSE  &
                                                                  is.na(labs_wl$int_fund) == FALSE &
                                                                  is.na(labs_wl$outgroup_pct)==F &
                                                                  is.na(labs_wl$natl)==F &
                                                                  is.na(labs_wl$total_domains)==F &
                                                                  is.na(labs_wl$DOMAIN_HUMA)==F &

                                                                  is.na(labs_wl$ORG_MISSION_INFO_2)==F]), df = 20)



full_reg<-lm(self_ind_pct~rebel+int_fund+natl+total_domains+DOMAIN_HUMA+ORG_MISSION_INFO_2+SEX+higher_edu+northern, data=labs_wl)
full_reg_self<-coeftest(full_reg, vcov =  sandwich::
                          vcovCL(full_reg, labs_wl$DEPT.x[is.na(labs_wl$rebel) == FALSE  &
                                                         is.na(labs_wl$int_fund) == FALSE &
                                                         is.na(labs_wl$self_ind_pct)==F &
                                                         is.na(labs_wl$natl)==F &
                                                         is.na(labs_wl$total_domains)==F &
                                                         is.na(labs_wl$DOMAIN_HUMA)==F &
                                                         is.na(labs_wl$SEX) == FALSE &
                                                         is.na(labs_wl$higher_edu)==F &
                                                         is.na(labs_wl$northern)==F &
                                                         is.na(labs_wl$ORG_MISSION_INFO_2)==F]), df = 20)

full_reg2<-lm(outgroup_pct~rebel+int_fund+natl+total_domains+DOMAIN_HUMA+ORG_MISSION_INFO_2+SEX+higher_edu+northern, data=labs_wl)
full_reg_out<-coeftest(full_reg2, vcov =  sandwich::
                         vcovCL(full_reg2, labs_wl$DEPT.x[is.na(labs_wl$rebel) == FALSE  &
                                                         is.na(labs_wl$int_fund) == FALSE &
                                                         is.na(labs_wl$outgroup_pct)==F &
                                                         is.na(labs_wl$natl)==F &
                                                         is.na(labs_wl$total_domains)==F &
                                                         is.na(labs_wl$DOMAIN_HUMA)==F &
                                                         is.na(labs_wl$SEX) == FALSE &
                                                         is.na(labs_wl$higher_edu)==F &
                                                         is.na(labs_wl$northern)==F &
                                                         is.na(labs_wl$ORG_MISSION_INFO_2)==F]), df = 20)





stargazer(base_self.lm2, self_ind.lm2, org_ctrls.self, full_reg_self,
          base_out.lm2, out.lm2, org_ctrls_out.out,full_reg_out, style = "apsr",
          type="text",
          column.labels   = c("Percent allocated to self", "Percent allocated to outgroups"),
          column.separate = c(4, 4),
          add.lines = list(c("Observations", paste(nobs(base_self.lm)),paste(nobs(self_ind.lm)),
                             paste(nobs(org_ctrls)),paste(nobs(full_reg2)),
                             paste(nobs(base_out.lm)),paste(nobs(out.lm)),
                             paste(nobs(org_ctrls_out)),paste(nobs(full_reg2))),
                           c("Individual controls", "No", "Yes","No", "Yes","No", "Yes", "No", "Yes"),
                           c("Organization controls", "No", "No","Yes", "Yes","No", "No", "Yes", "Yes")),
          keep = c("rebel"),
          covariate.labels = c("Wartime leader from rebel depts"),
          out = "TableA19.tex")
# Table A20: Results with SEs clustered at Session level ---------------------------------------------------------


rm(base_self.lm2, self_ind.lm2, org_ctrls.self, full_reg_self,
   base_out.lm2, out.lm2, org_ctrls_out.out,full_reg_out)



base_self.lm<-lm(self_ind_pct~rebel, data=labs)
self_ind.lm_s<-coeftest(base_self.lm, vcov =  sandwich::vcovCL(base_self.lm, labs$session[is.na(labs$rebel) == FALSE &
                                                                                            is.na(labs$self_ind_pct)==F]), df = 5)
base_out.lm<-lm(outgroup_pct~rebel, data=labs)
outgroup.lm_s<-coeftest(base_out.lm, vcov =  sandwich::vcovCL(base_out.lm, labs$session[is.na(labs$rebel) == FALSE]), df = 5)



stargazer(self_ind.lm_s, outgroup.lm_s,
          column.labels   = c("Percent allocated to self", "Percent allocated to outgroups"),
          covariate.labels = c("Civil society leader from rebel depts"),

          add.lines = list(c("Observations", paste(nobs(self_ind.lm_s)),paste(nobs(outgroup.lm_s)))),
          out = "TableA20.tex")



# Table A21: Outcomes for minority/northerners in Rebel Areas -------------------------------------------------------


TableA21<-data.frame(c("percent allocated to self", "percent allocated to outgroup"))
TableA21[2,2]<-t.test(labs$outgroup_pct[labs$minority==1 & labs$rebel==1],labs$outgroup_pct[labs$minority==0 & labs$rebel==1] )$estimate[1]
TableA21[2,3]<-t.test(labs$outgroup_pct[labs$minority==1 & labs$rebel==1],labs$outgroup_pct[labs$minority==0 & labs$rebel==1] )$estimate[2]
TableA21[2,4]<-t.test(labs$outgroup_pct[labs$minority==1 & labs$rebel==1],labs$outgroup_pct[labs$minority==0 & labs$rebel==1] )$p.value

TableA21[1,2]<-t.test(labs$self_ind_pct[labs$minority==1 & labs$rebel==1],labs$self_ind_pct[labs$minority==0 & labs$rebel==1] )$estimate[1]
TableA21[1,3]<-t.test(labs$self_ind_pct[labs$minority==1 & labs$rebel==1],labs$self_ind_pct[labs$minority==0 & labs$rebel==1] )$estimate[2]
TableA21[1,4]<-t.test(labs$self_ind_pct[labs$minority==1 & labs$rebel==1],labs$self_ind_pct[labs$minority==0 & labs$rebel==1] )$p.value
TableA21[1,5]<-sum(!is.na(labs$self_ind_pct))
TableA21[2,5]<-sum(!is.na(labs$outgroup_pct))

names(TableA21)<-c("outcome", "Minority group", "Northerner", "p-value", "N")
TableA21<- TableA21 %>% mutate_if(is.numeric, round, 2)

print(xtable(TableA21), include.rownames=F, file="TableA21.tex")

# Table A22: Outcomes for Christians/non-Christians in Rebel Areas -------------------------------------------------------


TableA22<-data.frame(c("percent allocated to self", "percent allocated to outgroup"))
TableA22[2,2]<-t.test(labs$outgroup_pct[labs$christian==1 & labs$rebel==1],labs$outgroup_pct[labs$christian==0 & labs$rebel==1] )$estimate[1]
TableA22[2,3]<-t.test(labs$outgroup_pct[labs$christian==1 & labs$rebel==1],labs$outgroup_pct[labs$christian==0 & labs$rebel==1] )$estimate[2]
TableA22[2,4]<-t.test(labs$outgroup_pct[labs$christian==1 & labs$rebel==1],labs$outgroup_pct[labs$christian==0 & labs$rebel==1] )$p.value

TableA22[1,2]<-t.test(labs$self_ind_pct[labs$christian==1 & labs$rebel==1],labs$self_ind_pct[labs$christian==0 & labs$rebel==1] )$estimate[1]
TableA22[1,3]<-t.test(labs$self_ind_pct[labs$christian==1 & labs$rebel==1],labs$self_ind_pct[labs$christian==0 & labs$rebel==1] )$estimate[2]
TableA22[1,4]<-t.test(labs$self_ind_pct[labs$christian==1 & labs$rebel==1],labs$self_ind_pct[labs$christian==0 & labs$rebel==1] )$p.value
TableA22[1,5]<-sum(!is.na(labs$self_ind_pct))
TableA22[2,5]<-sum(!is.na(labs$outgroup_pct))

names(TableA22)<-c("outcome", "Christians", "Non-Christians", "p-value", "N")
TableA22<- TableA22 %>% mutate_if(is.numeric, round, 2)

print(xtable(TableA22), include.rownames=F, file="TableA22.tex")


# Table A23: Post-treatment matching --------------------------------------------------------

labs_match<-subset(labs, !is.na(created_prewar) & !is.na(AGE) &
                      !is.na(EDUCATION) & !is.na(northern) &
                      !is.na(SEX) &!is.na(outgroup_pct) &!is.na(muslim) &!is.na(self_ind_pct))

Y<-labs_match$self_ind_pct
Y2<-labs_match$outgroup_pct
Tr<-labs_match$rebel

glm<-glm(Tr ~ AGE+SEX+EDUCATION+northern+created_prewar+muslim, family=binomial, data=labs_match)

rr_m_self<-Match(Y=Y, Tr=Tr, X=glm$fitted.values)
rr_m_outgroup<-Match(Y=Y2, Tr=Tr, X=glm$fitted.values)

TableA23<-data.frame(c("percent allocated to self","", "percent allocated to outgroup", ""))


TableA23[1,2]<-rr_m_self$est
TableA23[2,2]<-rr_m_self$se
TableA23[1,3]<-ifelse(((1 - pnorm(abs(rr_m_self$est/rr_m_self$se))) * 2)<.01,  "***",
                              ifelse(((1 - pnorm(abs(rr_m_self$est/rr_m_self$se))) * 2)<.05, "**",
                                     ifelse(((1 - pnorm(abs(rr_m_self$est/rr_m_self$se))) * 2)<.1, "*",NA)))


TableA23[3,2]<-rr_m_outgroup$est
TableA23[4,2]<-rr_m_outgroup$se
TableA23[3,3]<-ifelse(((1 - pnorm(abs(rr_m_outgroup$est/rr_m_outgroup$se))) * 2)<.01,  "***",
                              ifelse(((1 - pnorm(abs(rr_m_outgroup$est/rr_m_outgroup$se))) * 2)<.05, "**",
                                     ifelse(((1 - pnorm(abs(rr_m_outgroup$est/rr_m_outgroup$se))) * 2)<.1, "*",NA)))




names(TableA23)<-c("outcome", "ATT", "")
print(xtable(TableA23, digits=c(0,0,4,0)),include.rownames=F, file="TableA23.tex")



# Table A24: Analysis of population frame -------------------------------------------------

#this is the original sample frame, but is now anonymized to remove any names of orgs or leaders

#this is the individual level data
orig<-read.csv("original_sample_frame_anon.csv")

#this is aggregated at the org-level
orig2<-orig %>% distinct(ID2, .keep_all = TRUE)

orig2<-aggregate(orig2, list(orig2$ID2), mean)

TableA24<-data.frame(c("Northerner", "Woman", "Women's organization", "Health organization"))
TableA24[1,2]<-t.test(orig$northern[orig$rebel==1],orig$northern[orig$rebel==0])$estimate[1]
TableA24[1,3]<-t.test(orig$northern[orig$rebel==1],orig$northern[orig$rebel==0])$estimate[2]
TableA24[1,4]<-t.test(orig$northern[orig$rebel==1],orig$northern[orig$rebel==0])$p.value
TableA24[2,2]<-t.test(orig$female[orig$rebel==1],orig$female[orig$rebel==0])$estimate[1]
TableA24[2,3]<-t.test(orig$female[orig$rebel==1],orig$female[orig$rebel==0])$estimate[2]
TableA24[2,4]<-t.test(orig$female[orig$rebel==1],orig$female[orig$rebel==0])$p.value
TableA24[3,2]<-t.test(orig2$org_femme[orig2$rebel==1],orig2$org_femme[orig2$rebel==0])$estimate[1]
TableA24[3,3]<-t.test(orig2$org_femme[orig2$rebel==1],orig2$org_femme[orig2$rebel==0])$estimate[2]
TableA24[3,4]<-t.test(orig2$org_femme[orig2$rebel==1],orig2$org_femme[orig2$rebel==0])$p.value
TableA24[4,2]<-t.test(orig2$health[orig2$rebel==1],orig2$health[orig2$rebel==0])$estimate[1]
TableA24[4,3]<-t.test(orig2$health[orig2$rebel==1],orig2$health[orig2$rebel==0])$estimate[2]
TableA24[4,4]<-t.test(orig2$health[orig2$rebel==1],orig2$health[orig2$rebel==0])$p.value
TableA24[1,5]<-nrow(orig)
TableA24[2,5]<-nrow(orig)
TableA24[3,5]<-nrow(orig2)
TableA24[4,5]<-nrow(orig2)

names(TableA24)<-c("", "Rebel", "Govt", "p-value", "N")

print(xtable(TableA24), include.rownames=F, file="TableA24.tex")

# Table A25: Follow-up survey, summary stats --------------------------------------------------------

#df<-read.csv("Follow-up_survey_clean_12122022.csv")
df<-read.csv("Follow-up_survey.csv")

vars<-c("personal_org1", "personal_org2", "personal_org3","results_guess_ego_same",  "results_guess_disc_same",
        "results_guess_ego_correct","results_guess_disc_correct", "alt_surprised", "disc_surprised", "alt_agree2", "disc_agree2")

vars2<-c("personal only", "both perso/org", "org only", "Guess - altruism same","Guess - discrimination same ","Guess - altruism correct", "Guess - discrimination correct","Not surprised by Altruism finding", "Not surprised by Discrmination finding", "Altruism finding - agree", "Discrimination finding - agree")
TableA25<-t(stat.desc(df[vars]))

TableA25<-as.data.frame(TableA25[,c("mean", "min", "max", "std.dev", "nbr.val")])
TableA25$des<-vars2

TableA25<-TableA25[c(6,1,2,3,4,5)]

colnames(TableA25)<-c("description", "mean", "min", "max", "std.dev", "n")

print(xtable(TableA25, digits=c(0,0,2,0,0,2,0)), include.rownames=F, file="TableA25.tex")

# Table A26: Comparing followup Reb.vs govt -------------------------------

vars<-c("results_guess_ego_correct","results_guess_disc_correct", "alt_surprised", "disc_surprised")

vars2<-c("Guess - altruism correct", "Guess - discrimination correct","Not surprised by Altruism finding", "Not surprised by Discrmination finding")


TableA26<-data.frame((vars2))
for (i in 1:length(vars)){
  TableA26[i,2]<-lapply(df[vars[i]], function(x) {t.test(x[df$rebel==1], x[df$rebel==0])$estimate[1]})
  TableA26[i,3]<-lapply(df[vars[i]], function(x) {t.test(x[df$rebel==1], x[df$rebel==0])$estimate[2]})
  TableA26[i,4]<-lapply(df[vars[i]], function(x) {t.test(x[df$rebel==1], x[df$rebel==0])$p.value})
  TableA26[i,5]<-lapply(df[vars[i]], function(x) {min(x, na.rm=T)})
  TableA26[i,6]<-lapply(df[vars[i]], function(x) {max(x, na.rm=T)})
  TableA26[i,7]<-lapply(df[vars[i]], function(x) {sum(!is.na(x))})
}

names(TableA26)<-c("var", "Rebel","Govt", "p-value", "min", "max", "N")

print(xtable(TableA26), include.rownames=F, file="TableA26.tex")


# Table A27: Reasons for findings -----------------------------------------

reasons<-c("fear_reason", "poverty_reason", "aid_reason", "corrupt_reason", "disc_reason")
reasons2<-c("CSO leaders lived in uncertainty","CSO leaders are poorer", "CSO leaders get less aid", "CSO leaders are corrupt", "CSO leaders are victims of discrimination")

TableA27<-t(stat.desc(df[reasons]))

TableA27<-as.data.frame(TableA27[,c("mean", "min", "max", "std.dev", "nbr.val")])
TableA27$des<-reasons2

TableA27<-TableA27[c(6,1,2,3,4,5)]

colnames(TableA27)<-c("description", "mean", "min", "max", "std.dev", "n")

print(xtable(TableA27), digits=c(0,0,2,0,0,2,0),include.rownames=F, file="TableA27.tex")




